{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Graphical Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import fetch_openml\n",
    "from prml import bayesnet as bn\n",
    "\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = bn.discrete([0.1, 0.9])\n",
    "f = bn.discrete([0.1, 0.9])\n",
    "\n",
    "g = bn.discrete([[[0.9, 0.8], [0.8, 0.2]], [[0.1, 0.2], [0.2, 0.8]]], b, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.1 0.9])\n",
      "f: DiscreteVariable(proba=[0.1 0.9])\n",
      "g: DiscreteVariable(proba=[0.315 0.685])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "f: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(observed=[1. 0.])\n",
      "f: DiscreteVariable(proba=[0.11111111 0.88888889])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.3 Illustration: Image de-noising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/9s/lky4p_js2czgsr4_5962ffbw0000gn/T/ipykernel_11247/693879585.py:3: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.\n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  binarized_img = (x > 127).astype(np.int).reshape(28, 28)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x14909fc40>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAALOklEQVR4nO3dT6xmdX3H8fenqBsk6VDCZIpYbMPOBTaETUlDFxrKZnBhI6sxNrkuSmN3EruQxJiYprXLJmMkThuLMQEKIU2VECOuDAOhMDhRqBl1nMlMyLQRVyp8u7hnyHW4/+Y5z3nOc+f7fiVPnuc597nnfDncz/x+5/e75/5SVUi69v3e3AVIWg3DLjVh2KUmDLvUhGGXmnjPKg+WxKF/aWJVle22j2rZk9yb5EdJXk/y0Jh9SZpWFp1nT3Id8GPgo8BZ4Hnggar64S7fY8suTWyKlv0u4PWq+klV/Rr4JnB0xP4kTWhM2G8Bfr7l/dlh2+9IspHkZJKTI44laaQxA3TbdRXe1U2vquPAcbAbL81pTMt+Frh1y/sPAOfGlSNpKmPC/jxwe5IPJXkf8EngqeWUJWnZFu7GV9VvkzwIfBu4Dnikql5dWmWSlmrhqbeFDuY1uzS5SX6pRtLBYdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhOGXWrCsEtNGHapCcMuNWHYpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJwy41YdilJgy71IRhl5ow7FITCy/ZLI01dgXhZNvFSpe2/zHHXkejwp7kDPAm8Bbw26q6cxlFSVq+ZbTsf1FVbyxhP5Im5DW71MTYsBfwnSQvJNnY7gNJNpKcTHJy5LEkjZAxgxhJ/rCqziW5GXgG+Nuqem6Xz083YqIDxwG6aVTVtsWNatmr6tzwfBF4ArhrzP4kTWfhsCe5PskNl18DHwNOLaswScs1ZjT+MPDE0J15D/DvVfVfS6lKV2XK7uo66/rfvahR1+xXfTCv2SfhD/3qtbtml3RwGHapCcMuNWHYpSYMu9SEt7iugKPli1nnEe+DyJZdasKwS00YdqkJwy41YdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmvB+9hUYu3LJtbryiff5r5Ytu9SEYZeaMOxSE4ZdasKwS00YdqkJwy414Tz7Gpjz76N3PXZHe7bsSR5JcjHJqS3bbkzyTJLXhudD05Ypaaz9dOO/Dtx7xbaHgGer6nbg2eG9pDW2Z9ir6jng0hWbjwInhtcngPuXW5akZVv0mv1wVZ0HqKrzSW7e6YNJNoCNBY8jaUkmH6CrquPAcYAk3vkgzWTRqbcLSY4ADM8Xl1eSpCksGvangGPD62PAk8spR9JUso97qR8F7gFuAi4AXwD+A/gW8EHgZ8AnqurKQbzt9mU3fgIH9X52TaOqtv2ftmfYl8mwT8Owa6udwu6vy0pNGHapCcMuNWHYpSYMu9SEt7heA3YbMffPNesyW3apCcMuNWHYpSYMu9SEYZeaMOxSE4ZdasJ59mvc2OWex87Te9fc+rBll5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmnGdvbuw8/F52+37n4FfLll1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmnCeXbuach7ee+VXa8+WPckjSS4mObVl28NJfpHkpeFx37RlShprP934rwP3brP9n6vqjuHxn8stS9Ky7Rn2qnoOuLSCWiRNaMwA3YNJXh66+Yd2+lCSjSQnk5wccSxJI2U/gyRJbgOerqoPD+8PA28ABXwROFJVn97Hflxl8Boz58KRDtBtr6q2PTELtexVdaGq3qqqt4GvAneNKU7S9BYKe5IjW95+HDi102clrYc959mTPArcA9yU5CzwBeCeJHew2Y0/A3xmuhK1zsZ0pae8Vx7s5l9pX9fsSzuY1+zaYuqfva5hX+o1u6SDx7BLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhOGXWrCsEtNGHapCcMuNWHYpSb8U9IaZc6/VKOrY8suNWHYpSYMu9SEYZeaMOxSE4ZdasKwS004z96c8+R92LJLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhPOs1/jDvI8etdVWKeyZ8ue5NYk301yOsmrST47bL8xyTNJXhueD01frqRF7bk+e5IjwJGqejHJDcALwP3Ap4BLVfXlJA8Bh6rqc3vs6+A2MweULXs/C6/PXlXnq+rF4fWbwGngFuAocGL42Ak2/wGQtKau6po9yW3AR4AfAIer6jxs/oOQ5OYdvmcD2BhZp6SR9uzGv/PB5P3A94AvVdXjSf6vqn5/y9f/t6p2vW63G796duP7WbgbD5DkvcBjwDeq6vFh84Xhev7ydf3FZRQqaRr7GY0P8DXgdFV9ZcuXngKODa+PAU8uvzzBZuu86GNuSRZ+aLn2Mxp/N/B94BXg7WHz59m8bv8W8EHgZ8AnqurSHvua/6fvAFqH0C7K0K7eTt34fV+zL4NhX4xh19UYdc0u6eAz7FIThl1qwrBLTRh2qQlvcV2CgzxavhdH068dtuxSE4ZdasKwS00YdqkJwy41YdilJgy71ITz7INrea58N86j92HLLjVh2KUmDLvUhGGXmjDsUhOGXWrCsEtNtJlnv5bn0Z0r137YsktNGHapCcMuNWHYpSYMu9SEYZeaMOxSE/tZn/3WJN9NcjrJq0k+O2x/OMkvkrw0PO6bvtzFjVknfN0f0n7sZ332I8CRqnoxyQ3AC8D9wF8Bv6qqf9z3wVyyWZrcTks27/kbdFV1Hjg/vH4zyWngluWWJ2lqV3XNnuQ24CPAD4ZNDyZ5OckjSQ7t8D0bSU4mOTmuVElj7NmNf+eDyfuB7wFfqqrHkxwG3gAK+CKbXf1P77EPu/HSxHbqxu8r7EneCzwNfLuqvrLN128Dnq6qD++xH8MuTWynsO9nND7A14DTW4M+DNxd9nHg1NgiJU1nP6PxdwPfB14B3h42fx54ALiDzW78GeAzw2DebvuyZZcmNqobvyyGXZrewt14SdcGwy41YdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhOrXrL5DeCnW97fNGxbR+ta27rWBda2qGXW9kc7fWGl97O/6+DJyaq6c7YCdrGuta1rXWBti1pVbXbjpSYMu9TE3GE/PvPxd7Outa1rXWBti1pJbbNes0tanblbdkkrYtilJmYJe5J7k/woyetJHpqjhp0kOZPklWEZ6lnXpxvW0LuY5NSWbTcmeSbJa8PztmvszVTbWizjvcsy47Oeu7mXP1/5NXuS64AfAx8FzgLPAw9U1Q9XWsgOkpwB7qyq2X8BI8mfA78C/vXy0lpJ/gG4VFVfHv6hPFRVn1uT2h7mKpfxnqi2nZYZ/xQznrtlLn++iDla9ruA16vqJ1X1a+CbwNEZ6lh7VfUccOmKzUeBE8PrE2z+sKzcDrWthao6X1UvDq/fBC4vMz7rudulrpWYI+y3AD/f8v4s67XeewHfSfJCko25i9nG4cvLbA3PN89cz5X2XMZ7la5YZnxtzt0iy5+PNUfYt1uaZp3m//6sqv4U+Evgb4buqvbnX4A/YXMNwPPAP81ZzLDM+GPA31XVL+esZatt6lrJeZsj7GeBW7e8/wBwboY6tlVV54bni8ATbF52rJMLl1fQHZ4vzlzPO6rqQlW9VVVvA19lxnM3LDP+GPCNqnp82Dz7uduurlWdtznC/jxwe5IPJXkf8EngqRnqeJck1w8DJyS5HvgY67cU9VPAseH1MeDJGWv5HeuyjPdOy4wz87mbffnzqlr5A7iPzRH5/wH+fo4adqjrj4H/Hh6vzl0b8Cib3brfsNkj+mvgD4BngdeG5xvXqLZ/Y3Np75fZDNaRmWq7m81Lw5eBl4bHfXOfu13qWsl589dlpSb8DTqpCcMuNWHYpSYMu9SEYZeaMOxSE4ZdauL/AfayA45jFy/FAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, _ = fetch_openml(\"mnist_784\", return_X_y=True, as_frame=False)\n",
    "x = x[0]\n",
    "binarized_img = (x > 127).astype(np.int).reshape(28, 28)\n",
    "plt.imshow(binarized_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x14918b370>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAMpElEQVR4nO3dT8gc9R3H8c+n/rmo0KQ24WmMjS3ePGifh1wqxR6UNJfowaKniIXHQy32ptiDggihtJaeCrEG02IVwViDlGoQMZ4kT4KNiaEmlVRjHvJU0tJ4sppvDzuRx7h/nmdnZmdmv+8XLLs7z+7Od2f388xv5jc7P0eEAEy/rzVdAIDJIOxAEoQdSIKwA0kQdiCJSyc5M9tTuet/dnZ26N8PHjzY2PzrnndXNf2Z1Ski3G+6y3S92d4i6beSLpH0+4jYMeLxUxn2UcvQ7rvsJzL/uufdVU1/ZnWqPOy2L5H0nqRbJZ2SdEDS3RHx7pDnEPYJz7/LX9o6Nf2Z1WlQ2Mtss2+WdCIi3o+ITyU9J2lbidcDUKMyYd8g6cNl908V077E9rztBdsLJeYFoKQyO+j6NRW+0jaKiJ2SdkrT24wHuqDMmv2UpI3L7l8j6XS5cgDUpUzYD0i63vZ1ti+XdJekvdWUBaBqYzfjI+Iz2/dLekW9rrddEXG0sso6pOk9t8PmP817ncvI+L5L9bOvemZss08cYc+njq43AB1C2IEkCDuQBGEHkiDsQBKEHUhior9nb1LZLsY6u6jq7B6ja236DPu+zM3NDfwba3YgCcIOJEHYgSQIO5AEYQeSIOxAEmm63kZ1QTU5wCXdY1iNcb8vrNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IIk0/eyjTGtfd9M/7a1z0MkyPw3OeNZd1uxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kMTU9LPX3W9apr+4yT7dLvcXlz1GoMzzp7EfvlTYbZ+UdE7S55I+i4jBJ60G0Kgq1uw/jIiPK3gdADVimx1IomzYQ9Krtg/anu/3ANvzthdsL5ScF4ASXGYnhu1vRcRp2+sk7ZP0s4jYP+TxtZ3VkR107VTncmtSmz+ziOhbXKk1e0ScLq6XJL0oaXOZ1wNQn7HDbvsK21dduC3pNklHqioMQLXK7I1fL+nFojlzqaQ/RcRfyxTT5uZuy5ttTZcwlq7W3VWlttlXPbMR2+zTejKCJvcnoB4t/+df/TY7gO4g7EAShB1IgrADSRB2IIlW/cS1zB7ONu8dxXj4TKvFmh1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkmhVP3sZbf7V2yhdrh3dwZodSIKwA0kQdiAJwg4kQdiBJAg7kARhB5KYmn72NmPkk/7a/L6nEWt2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiiU/3sw/pl2/yb7zb/Xr3ueXf1M6tbE8tl5Jrd9i7bS7aPLJu21vY+28eL6zW1VAegMitpxj8tactF0x6S9FpEXC/pteI+gBYbGfaI2C/p7EWTt0naXdzeLen2assCULVxt9nXR8SiJEXEou11gx5oe17S/JjzAVCR2nfQRcROSTslyTa/fAAaMm7X2xnbM5JUXC9VVxKAOowb9r2Sthe3t0t6qZpyANTFK+gDflbSLZKulnRG0iOS/izpeUnXSvpA0p0RcfFOvH6vRTO+Zdp8DMAoZX4P3+b3VVZE9H1zI8NeJcLePoR9+gwKO4fLAkkQdiAJwg4kQdiBJAg7kAQ/cZ2AJvd4N3265jo/szZ/5m3Emh1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkuhUP3tXlR2yuUx/ct3DRddZe5d/kVenYctlbm5u4N9YswNJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEhM9u+zc3FwsLCwMLoY+2c5p8gyvbf7MGz6HAWeXBTIj7EAShB1IgrADSRB2IAnCDiRB2IEkGMW1Am3u761bk+eln+blWsbY/ey2d9lesn1k2bRHbX9k++3isrXKYgFUbyXN+Kclbekz/TcRcWNx+Uu1ZQGo2siwR8R+SWcnUAuAGpXZQXe/7cNFM3/NoAfZnre9YHvwQfEAareiHXS2N0l6OSJuKO6vl/SxpJD0mKSZiLh3Ba/DDropww669qn0hzARcSYiPo+I85KelLS5THEA6jdW2G3PLLt7h6Qjgx4LoB1Gnjfe9rOSbpF0te1Tkh6RdIvtG9Vrxp+UdF99JbZf5uZknecgQLU4qAaNKfvdy/xPdhhOXgEkR9iBJAg7kARhB5Ig7EASrRqyOfORaHWZ5mU6ze+tDqzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiCJVvWzd7VftM39vXXPm5+pjmfYcqvrM2PNDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJTDTss7OzioiBl1Hqem7Z59seeumyssttWpVdLk18X1izA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASjOI65brcF97kMQptPkfBKGOP4mp7o+3XbR+zfdT2A8X0tbb32T5eXK+pumgA1Rm5Zrc9I2kmIg7ZvkrSQUm3S7pH0tmI2GH7IUlrIuLBEa/V3dVMR7FmH0/KNXtELEbEoeL2OUnHJG2QtE3S7uJhu9X7BwCgpVZ1DjrbmyTdJOktSesjYlHq/UOwvW7Ac+YlzZesE0BJK95BZ/tKSW9Iejwi9tj+T0R8fdnf/x0RQ7fbacZPHs348aRsxkuS7cskvSDpmYjYU0w+U2zPX9iuX6qiUAD1GNmMd+9f2FOSjkXEE8v+tFfSdkk7iuuXaqlwmSZOv9v2ebddm9eAw3S17mFWsjf+ZklvSnpH0vli8sPqbbc/L+laSR9IujMizo54rVLf2rYGjrAPNo2habtBzfhOHVTT1sAR9sEI++SV2mYH0H2EHUiCsANJEHYgCcIOJNGqIZtHaeue3S7vLR+lrcscq8eaHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeS6FQ/e52mua98mKz96F0+E824WLMDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBJp+tmnuR99WJ/wNL/vMtrcj17XMQCs2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgiZFht73R9uu2j9k+avuBYvqjtj+y/XZx2Vp/ueOzPfRS9vlNXup835i8uj6zlYzPPiNpJiIO2b5K0kFJt0v6saRPIuJXq3gTrT3CI+PJDDCdBg3ZPPIIuohYlLRY3D5n+5ikDdWWB6Buq9pmt71J0k2S3iom3W/7sO1dttcMeM687QXbC+VKBVDGyGb8Fw+0r5T0hqTHI2KP7fWSPpYUkh5Tr6l/74jXoBkP1GxQM35FYbd9maSXJb0SEU/0+fsmSS9HxA0jXoewAzUbFPaV7I23pKckHVse9GLH3QV3SDpStkgA9VnJ3vibJb0p6R1J54vJD0u6W9KN6jXjT0q6r9iZN+y1OrtmH6bsWp9WRft0+TMp1YyvCmEfb95t/mJNqy5/JmM34wFMB8IOJEHYgSQIO5AEYQeSIOxAEhM9lfTs7KwWFuo5RL5sV0iTXSl1dt3RLTieaXxfrNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IIlJ/8T1X5L+uWzS1eqd2qqN2lpbW+uSqG1cVdb27Yj4Zr8/TDTsX5m5vRARc40VMERba2trXRK1jWtStdGMB5Ig7EASTYd9Z8PzH6attbW1LonaxjWR2hrdZgcwOU2v2QFMCGEHkmgk7La32P677RO2H2qihkFsn7T9TjEMdaPj0xVj6C3ZPrJs2lrb+2wfL677jrHXUG2tGMbbg4cZb3TZDalrIstt4tvsti+R9J6kWyWdknRA0t0R8e5ECxnA9klJcxHR+AEYtn8g6RNJf7gwtJbtX0o6GxE7in+UayLiwZbU9qhWOYx3TbUNGmb8HjW47Koc/nwcTazZN0s6ERHvR8Snkp6TtK2BOlovIvZLOnvR5G2Sdhe3d6v3ZZm4AbW1QkQsRsSh4vY5SReGGW902Q2payKaCPsGSR8uu39K7RrvPSS9avug7fmmi+lj/YVhtorrdQ3Xc7GRw3hP0kXDjLdm2Y0z/HlZTYS938m92tT/9/2I+J6kH0n6adFcxcr8TtJ31RsDcFHSr5ssphhm/AVJP4+I/zZZy3J96prIcmsi7KckbVx2/xpJpxuoo6+IOF1cL0l6Ub3NjjY5c2EE3eJ6qeF6vhARZyLi84g4L+lJNbjsimHGX5D0TETsKSY3vuz61TWp5dZE2A9Iut72dbYvl3SXpL0N1PEVtq8odpzI9hWSblP7hqLeK2l7cXu7pJcarOVL2jKM96BhxtXwsmt8+POImPhF0lb19sj/Q9IvmqhhQF3fkfS34nK06dokPates+5/6rWIfiLpG5Jek3S8uF7botr+qN7Q3ofVC9ZMQ7XdrN6m4WFJbxeXrU0vuyF1TWS5cbgskARH0AFJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEv8HnPtiZ5Z3e7oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "indices = np.random.choice(binarized_img.size, size=int(binarized_img.size * 0.1), replace=False)\n",
    "noisy_img = np.copy(binarized_img)\n",
    "noisy_img.ravel()[indices] = 1 - noisy_img.ravel()[indices]\n",
    "plt.imshow(noisy_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "markov_random_field = np.array([\n",
    "        [[bn.discrete([0.5, 0.5], name=f\"p(z_({i},{j}))\") for j in range(28)] for i in range(28)], \n",
    "        [[bn.DiscreteVariable(2) for _ in range(28)] for _ in range(28)]])\n",
    "a = 0.9\n",
    "b = 0.9\n",
    "pa = [[a, 1 - a], [1 - a, a]]\n",
    "pb = [[b, 1 - b], [1 - b, b]]\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    bn.discrete(pb, markov_random_field[0, i, j], out=markov_random_field[1, i, j], name=f\"p(x_({i},{j})|z_({i},{j}))\")\n",
    "    if i != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i + 1, j]], name=f\"p(z_({i},{j}), z_({i+1},{j}))\")\n",
    "    if j != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i, j + 1]], name=f\"p(z_({i},{j}), z_({i},{j+1}))\")\n",
    "    markov_random_field[1, i, j].observe(noisy_img[i, j], proprange=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x149605a30>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAALYklEQVR4nO3dT6xcZ3nH8e+vATYhUp2msdwQGlplxyI0UTaNqnQBSrNxWFCRlRGVLoumojsiuiASQkJVS5eVjIhwKxqElKSxoqoQRYiwQnGiNHGwIAG5YGzZjdyqYUVJni7ucXRx5v7xnJk5c+/z/UijmTl3PPP4+P78vud955w3VYWkg++3pi5A0moYdqkJwy41YdilJgy71MR7VvlhSRz6l5asqjJr+6iWPcl9SX6U5PUkD495L0nLlXnn2ZNcB/wY+ChwDngeeLCqfrjDn7Fll5ZsGS373cDrVfXTqvoV8E3g6Ij3k7REY8J+C/DzLc/PDdt+Q5KNJKeSnBrxWZJGGjNAN6ur8K5uelUdB46D3XhpSmNa9nPArVuefwA4P64cScsyJuzPA7cn+VCS9wGfBE4upixJizZ3N76qfp3kIeDbwHXAo1X16sIqk7RQc0+9zfVhHrNLS7eUL9VI2j8Mu9SEYZeaMOxSE4ZdasKwS02s9Hz2dbbbFGQyczZD2jds2aUmDLvUhGGXmjDsUhOGXWrCsEtNOPU2cGpNB50tu9SEYZeaMOxSE4ZdasKwS00YdqkJwy414Ty7DqydTlvu+L0KW3apCcMuNWHYpSYMu9SEYZeaMOxSE4ZdasJ5di3VKlcJvhZj69qP8/Sjwp7kLPAm8Bbw66q6axFFSVq8RbTsf1pVbyzgfSQtkcfsUhNjw17Ad5K8kGRj1guSbCQ5leTUyM+SNELGDFQk+b2qOp/kZuAZ4K+q6rkdXr+eozVamnUdoBtrnQfoqmpmcaNa9qo6P9xfAp4E7h7zfpKWZ+6wJ7k+yQ1XHgMfA04vqjBJizVmNP4w8OTQnXkP8C9V9e8LqeqAWfZy0Ae1q6zFGnXMfs0f1vSY3bAfPO2O2SXtH4ZdasKwS00YdqkJwy414SmuWlvrPOK9H9myS00YdqkJwy41YdilJgy71IRhl5ow7FITzrOvgWWfFSeBLbvUhmGXmjDsUhOGXWrCsEtNGHapCcMuNeE8+xpY53n0ZV75dp3/3geRLbvUhGGXmjDsUhOGXWrCsEtNGHapCcMuNeE8+wp0nk/u/HdfN7u27EkeTXIpyekt225M8kyS14b7Q8stU9JYe+nGfx2476ptDwPPVtXtwLPDc0lrbNewV9VzwOWrNh8FTgyPTwAPLLYsSYs27zH74aq6AFBVF5LcvN0Lk2wAG3N+jqQFWfoAXVUdB44DJNn5yoqSlmbeqbeLSY4ADPeXFleSpGWYN+wngWPD42PAU4spR9KyZA/XLH8MuBe4CbgIfAH4V+BbwAeBnwGfqKqrB/FmvZfd+BXzmvT9VNXMf9Rdw75Ihn31DHs/24Xdr8tKTRh2qQnDLjVh2KUmDLvUhKe4HnC7jbY7Wt+HLbvUhGGXmjDsUhOGXWrCsEtNGHapCcMuNeE8e3PLnkd3yeb1YcsuNWHYpSYMu9SEYZeaMOxSE4ZdasKwS004z65Rlnl1Ys+1XyxbdqkJwy41YdilJgy71IRhl5ow7FIThl1qwnl2TWaVKwhrDy17kkeTXEpyesu2R5L8IslLw+3+5ZYpaay9dOO/Dtw3Y/s/VNUdw+3fFluWpEXbNexV9RxweQW1SFqiMQN0DyV5eejmH9ruRUk2kpxKcmrEZ0kaKXsZJElyG/B0VX14eH4YeAMo4IvAkar69B7exxGZA2bKQTZPhJmtqmbumLla9qq6WFVvVdXbwFeBu8cUJ2n55gp7kiNbnn4cOL3dayWth13n2ZM8BtwL3JTkHPAF4N4kd7DZjT8LfGZ5JWqdjelKO8++Wns6Zl/Yh3nMri3G/u55zD7bQo/ZJe0/hl1qwrBLTRh2qQnDLjXhKa7at7zU9LWxZZeaMOxSE4ZdasKwS00YdqkJwy41YdilJpxn1yieprp/2LJLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhPOs2tHzqMfHLbsUhOGXWrCsEtNGHapCcMuNWHYpSYMu9SE8+wH3H6eJ/e674u1a8ue5NYk301yJsmrST47bL8xyTNJXhvuDy2/XEnz2nV99iRHgCNV9WKSG4AXgAeATwGXq+rLSR4GDlXV53Z5r/3bzOxTtuz9zL0+e1VdqKoXh8dvAmeAW4CjwInhZSfY/A9A0pq6pmP2JLcBHwF+AByuqguw+R9Ckpu3+TMbwMbIOiWNtGs3/p0XJu8Hvgd8qaqeSPI/VfXbW37+31W143G73fjVsxvfz9zdeIAk7wUeB75RVU8Mmy8Ox/NXjusvLaJQScuxl9H4AF8DzlTVV7b86CRwbHh8DHhq8eUJNlvneW9TSzL3TYu1l9H4e4DvA68Abw+bP8/mcfu3gA8CPwM+UVWXd3mv6X/79qF1CO28DO3qbdeN3/Mx+yIY9vkYdl2LUcfskvY/wy41YdilJgy71IRhl5rwFNc92s8j4mM4mn5w2LJLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhPOsw+cR9dBZ8suNWHYpSYMu9SEYZeaMOxSE4ZdasKwS004z37AOY+uK2zZpSYMu9SEYZeaMOxSE4ZdasKwS00YdqmJvazPfmuS7yY5k+TVJJ8dtj+S5BdJXhpu9y+/3PU0Zg3yZd+kK/ayPvsR4EhVvZjkBuAF4AHgz4FfVtXf7fnD1njJ5jEXrzBUWifbLdm86zfoquoCcGF4/GaSM8Atiy1P0rJd0zF7ktuAjwA/GDY9lOTlJI8mObTNn9lIcirJqXGlShpj1278Oy9M3g98D/hSVT2R5DDwBlDAF9ns6n96l/ewGy8t2Xbd+D2FPcl7gaeBb1fVV2b8/Dbg6ar68C7vY9ilJdsu7HsZjQ/wNeDM1qAPA3dXfBw4PbZIScuzl9H4e4DvA68Abw+bPw88CNzBZjf+LPCZYTBvp/faty27rbf2i1Hd+EUx7NLyzd2Nl3QwGHapCcMuNWHYpSYMu9SEYZeaWGnY77zzTqpqKbexPFVUB50tu9SEYZeaMOxSE4ZdasKwS00YdqkJwy41sepTXP8L+M8tm25i89JW62hda1vXusDa5rXI2n6/qn531g9WGvZ3fXhyqqrumqyAHaxrbetaF1jbvFZVm914qQnDLjUxddiPT/z5O1nX2ta1LrC2ea2ktkmP2SWtztQtu6QVMexSE5OEPcl9SX6U5PUkD09Rw3aSnE3yyrAM9aTr0w1r6F1KcnrLthuTPJPkteF+5hp7E9W2Fst477DM+KT7burlz1d+zJ7kOuDHwEeBc8DzwINV9cOVFrKNJGeBu6pq8i9gJPkT4JfAP11ZWivJ3wKXq+rLw3+Uh6rqc2tS2yNc4zLeS6ptu2XGP8WE+26Ry5/PY4qW/W7g9ar6aVX9CvgmcHSCOtZeVT0HXL5q81HgxPD4BJu/LCu3TW1roaouVNWLw+M3gSvLjE+673aoayWmCPstwM+3PD/Heq33XsB3kryQZGPqYmY4fGWZreH+5onrudquy3iv0lXLjK/Nvptn+fOxpgj7rAu6rdP83x9X1R8Bfwb85dBd1d78I/CHbK4BeAH4+ymLGZYZfxz466r63ylr2WpGXSvZb1OE/Rxw65bnHwDOT1DHTFV1fri/BDzJ5mHHOrl4ZQXd4f7SxPW8o6ouVtVbVfU28FUm3HfDMuOPA9+oqieGzZPvu1l1rWq/TRH254Hbk3woyfuATwInJ6jjXZJcPwyckOR64GOs31LUJ4Fjw+NjwFMT1vIb1mUZ7+2WGWfifTf58ufLurTzLpd9vp/NEfmfAH8zRQ3b1PUHwH8Mt1enrg14jM1u3f+x2SP6C+B3gGeB14b7G9eotn9mc2nvl9kM1pGJaruHzUPDl4GXhtv9U++7HepayX7z67JSE36DTmrCsEtNGHapCcMuNWHYpSYMu9SEYZea+H8p/FHxlv9bMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(10000):\n",
    "    i, j = np.random.choice(28, 2)\n",
    "    markov_random_field[1, i, j].send_message(proprange=3)\n",
    "restored_img = np.zeros_like(noisy_img)\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    restored_img[i, j] = np.argmax(markov_random_field[0, i, j].proba)\n",
    "plt.imshow(restored_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
